This question involves the use of simple linear regression on the Auto data set.
library(data.table)
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.5 ✓ purrr 0.3.4
## ✓ tibble 3.1.6 ✓ dplyr 1.0.7
## ✓ tidyr 1.1.4 ✓ stringr 1.4.0
## ✓ readr 2.1.1 ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::between() masks data.table::between()
## x dplyr::filter() masks stats::filter()
## x dplyr::first() masks data.table::first()
## x dplyr::lag() masks stats::lag()
## x dplyr::last() masks data.table::last()
## x purrr::transpose() masks data.table::transpose()
auto <- fread("../Assignment 1/Auto.csv")
auto$horsepower <- as.numeric(auto$horsepower)
## Warning: NAs introduced by coercion
Use the
lm()function to perform a simple linear regression withmpgas the response andhorsepoweras the predictor. Use thesummary()function to print the results. Comment on the output. For example:
lm <- lm(mpg ~ horsepower, data = auto)
lm
##
## Call:
## lm(formula = mpg ~ horsepower, data = auto)
##
## Coefficients:
## (Intercept) horsepower
## 39.9359 -0.1578
summary(lm)
##
## Call:
## lm(formula = mpg ~ horsepower, data = auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.5710 -3.2592 -0.3435 2.7630 16.9240
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 39.935861 0.717499 55.66 <2e-16 ***
## horsepower -0.157845 0.006446 -24.49 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.906 on 390 degrees of freedom
## (5 observations deleted due to missingness)
## Multiple R-squared: 0.6059, Adjusted R-squared: 0.6049
## F-statistic: 599.7 on 1 and 390 DF, p-value: < 2.2e-16
Is there a relationship between the predictor and the re- sponse?
Yes. After horsepower is converted to numeric, there is a relationship between horsepower and mpg.
How strong is the relationship between the predictor and the response?
Very strong. The p-value is tiny.
Is the relationship between the predictor and the response positive or negative?
Negative. As horsepower increases, mpg decreases.
What is the predicted
mpgassociated with ahorsepowerof 98? What are the associated 95% confidence and prediction intervals?
predict(lm, auto)[98]
## 98
## 23.36216
# confidence interval
predict(lm, auto, interval = "confidence")[98,]
## fit lwr upr
## 23.36216 22.87497 23.84936
# prediction interval
predict(lm, auto, interval = "prediction")[98,]
## fit lwr upr
## 23.36216 13.70483 33.01950
Plot the response and the predictor. Use the
abline()function to display the least squares regression line.
plot(auto$horsepower, auto$mpg)
abline(lm)
Use the
plot()function to produce diagnostic plots of the least squares regression fit. Comment on any problems you see with the fit.
plot(lm)
From the textbook, Figure 3.9: “Plots of residuals versus predicted (or fitted) values for the Auto data set. In each plot, the red line is a smooth fit to the residuals, intended to make it easier to identify a trend. Left: A linear regression of mpg on horsepower. A strong pattern in the residuals indicates non-linearity in the data. Right: A linear regression of mpg on horsepower and horsepower. There is little pattern in the residuals.”
This question involves the use of multiple linear regression on the
Autodata set.
Produce a scatterplot matrix which includes all of the variables in the data set.
auto$name <- as.factor(auto$name)
pairs(auto)
Compute the matrix of correlations between the variables using the function
cor(). You will need to exclude thenamevariable, which is qualitative.
cor(auto[, -9])
## mpg cylinders displacement horsepower weight
## mpg 1.0000000 -0.7762599 -0.8044430 NA -0.8317389
## cylinders -0.7762599 1.0000000 0.9509199 NA 0.8970169
## displacement -0.8044430 0.9509199 1.0000000 NA 0.9331044
## horsepower NA NA NA 1 NA
## weight -0.8317389 0.8970169 0.9331044 NA 1.0000000
## acceleration 0.4222974 -0.5040606 -0.5441618 NA -0.4195023
## year 0.5814695 -0.3467172 -0.3698041 NA -0.3079004
## origin 0.5636979 -0.5649716 -0.6106643 NA -0.5812652
## acceleration year origin
## mpg 0.4222974 0.5814695 0.5636979
## cylinders -0.5040606 -0.3467172 -0.5649716
## displacement -0.5441618 -0.3698041 -0.6106643
## horsepower NA NA NA
## weight -0.4195023 -0.3079004 -0.5812652
## acceleration 1.0000000 0.2829009 0.2100836
## year 0.2829009 1.0000000 0.1843141
## origin 0.2100836 0.1843141 1.0000000
Use the
lm()function to perform a multiple linear regression withmpgas the response and all other variables exceptnameas the predictors. Use thesummary()function to print the results. Comment on the output. For instance:
lm <- lm(mpg ~ cylinders + displacement + horsepower + weight + acceleration + year + origin, auto)
lm
##
## Call:
## lm(formula = mpg ~ cylinders + displacement + horsepower + weight +
## acceleration + year + origin, data = auto)
##
## Coefficients:
## (Intercept) cylinders displacement horsepower weight
## -17.218435 -0.493376 0.019896 -0.016951 -0.006474
## acceleration year origin
## 0.080576 0.750773 1.426140
summary(lm)
##
## Call:
## lm(formula = mpg ~ cylinders + displacement + horsepower + weight +
## acceleration + year + origin, data = auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.5903 -2.1565 -0.1169 1.8690 13.0604
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -17.218435 4.644294 -3.707 0.00024 ***
## cylinders -0.493376 0.323282 -1.526 0.12780
## displacement 0.019896 0.007515 2.647 0.00844 **
## horsepower -0.016951 0.013787 -1.230 0.21963
## weight -0.006474 0.000652 -9.929 < 2e-16 ***
## acceleration 0.080576 0.098845 0.815 0.41548
## year 0.750773 0.050973 14.729 < 2e-16 ***
## origin 1.426141 0.278136 5.127 4.67e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.328 on 384 degrees of freedom
## (5 observations deleted due to missingness)
## Multiple R-squared: 0.8215, Adjusted R-squared: 0.8182
## F-statistic: 252.4 on 7 and 384 DF, p-value: < 2.2e-16
Is there a relationship between the predictors and the re- sponse?
There is a relationship between mpg and displacement, weight, year, and origin.
Which predictors appear to have a statistically significant relationship to the response?
There is a relationship between mpg and displacement, weight, year, and origin.
What does the coefficient for the
yearvariable suggest?
This suggests that the mpg increases alongside the year. (Cars become more efficient over time).
Use the
plot()function to produce diagnostic plots of the linear regression fit. Comment on any problems you see with the fit. Do the residual plots suggest any unusually large outliers? Does the leverage plot identify any observations with unusually high leverage?
plot(lm)
There are a decent amount of outliers among the dataset. There are a few observations with strangely high leverage.
Use the
*and:symbols to fit linear regression models with interaction effects. Do any interactions appear to be statistically significant?
summary(lm(mpg ~ cylinders * displacement + horsepower:weight, auto))
##
## Call:
## lm(formula = mpg ~ cylinders * displacement + horsepower:weight,
## data = auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.8087 -2.2841 -0.5424 2.1259 17.5877
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.524e+01 2.384e+00 23.167 < 2e-16 ***
## cylinders -3.831e+00 5.339e-01 -7.176 3.68e-12 ***
## displacement -1.388e-01 1.511e-02 -9.184 < 2e-16 ***
## cylinders:displacement 1.923e-02 2.171e-03 8.857 < 2e-16 ***
## horsepower:weight -2.231e-05 2.961e-06 -7.533 3.54e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.164 on 387 degrees of freedom
## (5 observations deleted due to missingness)
## Multiple R-squared: 0.7182, Adjusted R-squared: 0.7153
## F-statistic: 246.6 on 4 and 387 DF, p-value: < 2.2e-16
I guessed and it turns out that my suspicions were correct. cylinders:diplacement and horsepower:weight were significant.
Try a few different transformations of the variables, such as \(log(X)\), \(\sqrt{X}\), \(X^2\). Comment on your findings.
summary(lm(log(mpg) ~ cylinders^2, auto))
##
## Call:
## lm(formula = log(mpg) ~ cylinders^2, data = auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.61694 -0.12149 -0.00995 0.12358 0.62573
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.002754 0.032379 123.62 <2e-16 ***
## cylinders -0.165149 0.005664 -29.16 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1918 on 395 degrees of freedom
## Multiple R-squared: 0.6828, Adjusted R-squared: 0.682
## F-statistic: 850.3 on 1 and 395 DF, p-value: < 2.2e-16
summary(lm(log(mpg) ~ cylinders^2 + log(horsepower)^2:sqrt(weight), auto))
##
## Call:
## lm(formula = log(mpg) ~ cylinders^2 + log(horsepower)^2:sqrt(weight),
## data = auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.47818 -0.10004 -0.00615 0.10030 0.51869
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.4881118 0.0421172 106.563 <2e-16 ***
## cylinders -0.0211598 0.0107477 -1.969 0.0497 *
## log(horsepower):sqrt(weight) -0.0050924 0.0003451 -14.756 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1535 on 389 degrees of freedom
## (5 observations deleted due to missingness)
## Multiple R-squared: 0.7972, Adjusted R-squared: 0.7961
## F-statistic: 764.4 on 2 and 389 DF, p-value: < 2.2e-16
I’m a little surprised that the regressions didn’t break as easily as I thought they would with some basic transformations.
This question should be answered using the
Carseatsdata set.
carseats <- ISLR2::Carseats
Fit a multiple regression model to predict
SalesusingPrice,Urban, andUS.
lm <- lm(Sales ~ Price + Urban + US, carseats)
lm
##
## Call:
## lm(formula = Sales ~ Price + Urban + US, data = carseats)
##
## Coefficients:
## (Intercept) Price UrbanYes USYes
## 13.04347 -0.05446 -0.02192 1.20057
summary(lm)
##
## Call:
## lm(formula = Sales ~ Price + Urban + US, data = carseats)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.9206 -1.6220 -0.0564 1.5786 7.0581
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.043469 0.651012 20.036 < 2e-16 ***
## Price -0.054459 0.005242 -10.389 < 2e-16 ***
## UrbanYes -0.021916 0.271650 -0.081 0.936
## USYes 1.200573 0.259042 4.635 4.86e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.472 on 396 degrees of freedom
## Multiple R-squared: 0.2393, Adjusted R-squared: 0.2335
## F-statistic: 41.52 on 3 and 396 DF, p-value: < 2.2e-16
Provide an interpretation of each coefficient in the model. Be careful—some of the variables in the model are qualitative!
As Price increases by 1, Sales decrease by \(.054459\) When Urban is Yes, Sales decrease by \(.021916\) (but this isn’t significant) When US is Yes, Sales increases by \(1.021916\)
Write out the model in equation form, being careful to handle the qualitative variables properly.
\(\textrm{Sales}=13.04-.05\textrm{Price}-(.02\textrm{Urban}_{\textrm{Yes}})+(1.20\textrm{US}_{\textrm{Yes}})\)
For which of the predictors can you reject the null hypothesis \(H0:β_j=0\) ?
Price and USYes are significant enough with a p-value less than 0.05 so we are 95% confident that they are significant.
On the basis of your response to the previous question, fit a smaller model that only uses the predictors for which there is evidence of association with the outcome.
lm <- lm(Sales ~ Price + US, carseats)
lm
##
## Call:
## lm(formula = Sales ~ Price + US, data = carseats)
##
## Coefficients:
## (Intercept) Price USYes
## 13.03079 -0.05448 1.19964
summary(lm)
##
## Call:
## lm(formula = Sales ~ Price + US, data = carseats)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.9269 -1.6286 -0.0574 1.5766 7.0515
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.03079 0.63098 20.652 < 2e-16 ***
## Price -0.05448 0.00523 -10.416 < 2e-16 ***
## USYes 1.19964 0.25846 4.641 4.71e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.469 on 397 degrees of freedom
## Multiple R-squared: 0.2393, Adjusted R-squared: 0.2354
## F-statistic: 62.43 on 2 and 397 DF, p-value: < 2.2e-16
How well do the models in (a) and (e) fit the data?
The model in (e) fits much better than the one in (a) because it only has significant variables.
Using the model from (e), obtain \(95\%\) confidence intervals for the coefficient(s).
confint(lm)
## 2.5 % 97.5 %
## (Intercept) 11.79032020 14.27126531
## Price -0.06475984 -0.04419543
## USYes 0.69151957 1.70776632
Is there evidence of outliers or high leverage observations in the model from (e)?
par(mfrow = c(2,2))
plot(lm)
There are a few points that look to be outside the normal range. Points 51, 69, and 377 seem to be outliers and points 26, 50, and 368 seem to have high leverage. I’m not worried about them because there are only three in a good-sized dataset.
From the data available at http://www.stat.ufl.edu/~winner/datasets.html, obtain the data on “Fibre Diameters and Breaking Strenghs for Nextel 610 Fibres.” (please note that there is a typo on the website. It should be Strength and not Strengh) According to the description available there, the expectation is that the log of breaking strength of the fibre should be negatively and linearly related to diameter. (Note log here means natural log if not specified.)
dt <- fread("fibre.csv")
colnames(dt) <- c("stren", "diam")
head(dt)
## stren diam
## 1: 2944 11.0
## 2: 3206 11.3
## 3: 2751 12.0
## 4: 3089 11.3
## 5: 3406 11.0
## 6: 3275 11.2
Produce a scatter plot of breaking strength against diameter.
plot(dt$diam, dt$stren)
Produce a scatter plot of the log of breaking strength against diameter.
dt$logStren <- log(dt$stren)
plot(dt$logStren, dt$diam)
Produce a scatter plot of the log of breaking strength against the log of diameter.
dt$logDiam <- log(dt$diam)
plot(dt$logStren, dt$logDiam)
Regress breaking strength on diameter.
lm <- lm(stren ~ diam, data = dt)
lm
##
## Call:
## lm(formula = stren ~ diam, data = dt)
##
## Coefficients:
## (Intercept) diam
## 5312.8 -197.7
summary(lm)
##
## Call:
## lm(formula = stren ~ diam, data = dt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -603.46 -246.10 -6.76 223.64 1047.22
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5312.78 1040.41 5.106 5.61e-06 ***
## diam -197.73 91.91 -2.151 0.0365 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 335.7 on 48 degrees of freedom
## Multiple R-squared: 0.08794, Adjusted R-squared: 0.06894
## F-statistic: 4.628 on 1 and 48 DF, p-value: 0.03651
Regress the log of breaking strength on diameter.
lm <- lm(logStren ~ diam, data = dt)
lm
##
## Call:
## lm(formula = logStren ~ diam, data = dt)
##
## Coefficients:
## (Intercept) diam
## 8.76091 -0.06504
summary(lm)
##
## Call:
## lm(formula = logStren ~ diam, data = dt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.211992 -0.075119 0.003561 0.076367 0.293763
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.76091 0.33575 26.094 <2e-16 ***
## diam -0.06504 0.02966 -2.193 0.0332 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1083 on 48 degrees of freedom
## Multiple R-squared: 0.09105, Adjusted R-squared: 0.07211
## F-statistic: 4.808 on 1 and 48 DF, p-value: 0.0332
Regress the log of breaking strength on the log of diameter.
lm <- lm(logStren ~ logDiam, data = dt)
lm
##
## Call:
## lm(formula = logStren ~ logDiam, data = dt)
##
## Coefficients:
## (Intercept) logDiam
## 9.8328 -0.7454
summary(lm)
##
## Call:
## lm(formula = logStren ~ logDiam, data = dt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.21123 -0.07610 0.00396 0.07673 0.29398
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.8328 0.8209 11.978 4.99e-16 ***
## logDiam -0.7454 0.3385 -2.202 0.0325 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1083 on 48 degrees of freedom
## Multiple R-squared: 0.09174, Adjusted R-squared: 0.07282
## F-statistic: 4.848 on 1 and 48 DF, p-value: 0.03251
From the data available at http://www.stat.ufl.edu/~winner/datasets.html, obtain the data on “Variables associated with Permeability and Porosity of Rocks”
dt <- fread("rocks.csv")
colnames(dt) <- c("type", "density", "porosity", "logPermeability", "residue", "carbonate", "aGrain", "aSD", "bGrain", "bSD", "calcite", "dolomite")
head(dt)
## type density porosity logPermeability residue carbonate aGrain aSD bGrain
## 1: 1 2.6975 0.83 -1.2218 4.43 95.57 11.88 1.563 12.19
## 2: 1 2.6415 0.59 -1.6989 7.19 92.81 10.31 1.410 10.85
## 3: 1 2.6457 0.67 -1.6989 2.17 97.83 12.51 1.690 12.92
## 4: 1 2.6342 0.50 -1.6989 6.34 93.66 10.76 1.327 11.40
## 5: 1 2.6236 1.10 -1.2218 2.24 97.76 8.62 1.560 8.92
## 6: 2 2.7278 1.10 -1.6989 8.85 91.16 9.85 0.848 10.25
## bSD calcite dolomite
## 1: 1.379 74.44 21.08
## 2: 1.361 85.80 7.00
## 3: 1.759 93.71 4.13
## 4: 1.361 95.55 2.17
## 5: 1.427 77.38 20.41
## 6: 0.815 1.22 89.93
Fit a multiple regression model to predict porosity. Please provide a clean model with only significant variables. Interpret your model results and diagnostics.
lm <- lm(porosity ~ ., data = dt)
lmSumm <- summary(lm)
coefs <- coef(lmSumm)[-1,4] <= .05
# This prunes the regression until only significant variables remain
while(length(coefs) != sum(coefs)) {
form <- paste("porosity ~", paste(names(coefs[coefs]), "+", collapse = " "), collapse = " ")
form <- substr(form, 1, nchar(form) - 2)
form
lm <- lm(formula = form, data = dt)
lmSumm <- summary(lm)
coefs <- coef(lmSumm)[-1,4] <= .05
}
lm
##
## Call:
## lm(formula = form, data = dt)
##
## Coefficients:
## (Intercept) logPermeability
## 6.240 2.666
lmSumm
##
## Call:
## lm(formula = form, data = dt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.1522 -1.0096 -0.2271 0.3610 5.8290
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.240 2.046 3.050 0.00485 **
## logPermeability 2.666 1.222 2.181 0.03740 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.689 on 29 degrees of freedom
## Multiple R-squared: 0.141, Adjusted R-squared: 0.1113
## F-statistic: 4.759 on 1 and 29 DF, p-value: 0.0374
par(mfrow = c(2,2))
plot(lm)
This is a little silly, but the only significant variable is logPermeability and the intercept is even more significant than that is. The \(R^2\) value is far from great. I think that because the dataset is rather small, there are more issues that are more prevalent. This is seen in the diagnostic charts where the values and lines are all over the place. Points 5, 8, 10, and 20 seem to be the biggest problems, but that’s 13% of the total dataset.
Fit a multiple regression model to predict log(permeability). Please provide a clean model with only significant variables. Interpret your model results and diagnostics.
lm <- lm(logPermeability ~ ., data = dt)
lmSumm <- summary(lm)
coefs <- coef(lmSumm)[-1,4] <= .05
# This prunes the regression until only significant variables remain
while(length(coefs) != sum(coefs)) {
form <- paste("logPermeability ~", paste(names(coefs[coefs]), "+", collapse = " "), collapse = " ")
form <- substr(form, 1, nchar(form) - 2)
form
lm <- lm(formula = form, data = dt)
lmSumm <- summary(lm)
coefs <- coef(lmSumm)[-1,4] <= .05
}
lm
##
## Call:
## lm(formula = form, data = dt)
##
## Coefficients:
## (Intercept) porosity aGrain bGrain calcite
## -0.448085 0.063658 0.897464 -0.994242 0.003091
lmSumm
##
## Call:
## lm(formula = form, data = dt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.58844 -0.06918 0.02635 0.09975 0.40127
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.448085 0.499277 -0.897 0.37771
## porosity 0.063658 0.023650 2.692 0.01227 *
## aGrain 0.897464 0.430584 2.084 0.04710 *
## bGrain -0.994242 0.431684 -2.303 0.02952 *
## calcite 0.003091 0.001066 2.899 0.00751 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2019 on 26 degrees of freedom
## Multiple R-squared: 0.4452, Adjusted R-squared: 0.3598
## F-statistic: 5.216 on 4 and 26 DF, p-value: 0.003204
par(mfrow = c(2,2))
plot(lm)
This regression is much better, but still not particularly great. The \(R^2\) is much better and the Intercept is no longer significant. Porosity is significant which is expected based on the results of 5a. There are again a few outliers and leverage points which skew the results a bit.